%%% Written by  Daniel Lewis, 2020, for "Robust Inference in Models 
%%% Identified via Heteroskedasticity", Review of Economics and Statistics.
%
% This function computes robust confidence sets using the subset inference
% approach based on Theorem 3. It inverts the test statistic over a grid of
% values for the parameter of interest. It takes as arguments the estimated
% parameter vector, the outer product of residuals in each regime, an
% index for the parameter of interest, a grid over which to perturb the
% parameter of interest, and the level of the test. It outputs a robust
% confidence interval for the parameter of interest, and optionally
% produces confidence intervals under a set of distortions for an Andrews 
% (2018) identification test and the full grid of test statistics.
%
% Inputs:
% theta:    parameter vector
% etasq1:   outer product of regime 1 residuals (T1x3 unique vectorization)
% etasq2:   outer product of regime 2 residuals (T2x3 unique vectorization)
% pind:     index for parameter of interest
% range:    relative grid for parameter of interest
% level:    level of the test
% Outputs:
% CIbnds:   bounds of confidence set of level "level"
% DistCI:   bounds of confidence set for range of distortions rel. to level
% stats:    full grid of test statistics

function [CIbnds,DistCI,stats]=InvertTest(theta,etasq1,etasq2,pind,range,level)
%% Setup
estsub=theta;
estsub(pind)=[]; % nuisance parameters (alpha in the text)
theta0=theta(pind); % estimate for parameter of interst.
start1=estsub; % initiailze start values for search to right
% start1(5)=start1(5)-start1(3)*start1(4)/start1(2);
start2=estsub; % initiailze start values for search to left
% start2(5)=start2(5)-start2(3)*start2(4)/start2(2);
options=optimoptions(@fmincon,'Display','off'); % suppress output
warning('off','all')
Hrange=range+theta0; % center grid at estimate
b1=NaN(2,1); % will collect test statistics to right
b2=NaN(2,1); % will collect test statistics to left
lg=(length(Hrange)-1)/2; % length of grid in each direct
interval=ceil(lg/10); % interval at which to print output
T=length(etasq1)+length(etasq2); % Number of observations

%% Compute test statistics
for i=1:lg % looping over grid
    Htest=Hrange(lg+1+i); % select next grid value to the right
    f=@(thetahat) momfastsub(etasq1,etasq2,T,thetahat,Htest,pind); % set function for minimization
    try
        [est1,val,flag]=fmincon(f,start1,[],[],[],[],[-Inf eps*ones(1,4)],[],@ordering,options); % minimize, using previous argmin as start value, save test statistic in b1
        if flag>=0 && val<=f(start1) % if valid exit flag, save value of test statistic
            b1(i)=val;
            start1=est1; % save new argmin as start value for move to the right
        else
            b1(i)=f(start1);
        end
    catch
        b1(i)=f(start1); % if some numerical/convergence issue, save value at previous optimum
    end
    Htest=Hrange(lg+1-i); % select next grid value to the left
    f=@(thetahat) momfastsub(etasq1,etasq2,T,thetahat,Htest,pind); % set function for minimization
    try
        [est2,val,flag]=fmincon(f,start2,[],[],[],[],[-Inf eps*ones(1,4)],[],@ordering,options); % minimize, using previous argmin as start value, save test statistic in b2
        if flag>=0 && val<=f(start2) % if valid exit flag, save value of test statistic
            b2(i)=val;
            start2=est2; % save new argmin as start value for move to the left
        else
            b2(i)=f(start2);
        end
    catch
        b2(i)=f(start2); % if some numerical/convergence issue, save value at previous optimum
    end
    
    if interval*floor(i/interval)==i % ever gridlength/10 iterations, print summaries of the past 50 test statistics
        fprintf('Grid value = %d\n',i);
        mnr=min(b1(end-interval+1:end)); mxr=max(b1(end-interval+1:end)); % min and max test statistics to right
        mnl=min(b2(end-interval+1:end)); mxl=max(b2(end-interval+1:end)); % min and max test statistics to left
        fprintf('Minimum and maximum test statistics to the right: = %.2f, %.2f\n',mnr,mxr);
        fprintf('Minimum and maximum test statistics to the left: = %.2f, %.2f\n',mnl,mxl);
    end
    
end

stats=[flipud(b2);0;b1]; % combine test statistics in vector in order of grid values (0 for the actual estimate since model just-identified)
for j=2:length(stats)-1 % now if there are any NaNs in the test statistics, take midpoint of adjoining values
    if isnan(stats(j))
        stats(j)=(stats(j-1)+stats(j+1))/2;
    end
end
%% Compute confidence interval
lowind=find(b2<chi2inv(1-level,1)); % find grid values to the left the test does not reject
highind=find(b1<chi2inv(1-level,1)); % find grid values to the right the test does not reject

% In case there are no values found that have test statistics below the
% critical value:
if isempty(highind)
    highind=0;
end
if isempty(lowind)
    lowind=0;
end
CIbnds=[Hrange(lg-max(lowind)+1)  Hrange(max(highind)+lg+1)]; % find the last points in each direction for which the test does not reject

%% Compute intervals for given distortions for Andrews identification test
if nargout>1
    distortions=[.05 .1 .15 .2]; % Set of size distortion thresholds to evaluate
    dlevels=level+distortions; % computing corresponding levels
    DistCI=NaN(length(dlevels),2);
    for i=1:length(dlevels)
        lowind=find(b2<chi2inv(1-dlevels(i),1)); % find grid values to the left the test does not reject
        highind=find(b1<chi2inv(1-dlevels(i),1)); % find grid values to the right the test does not reject
        if isempty(highind)
            highind=0;
        end
        if isempty(lowind)
            lowind=0;
        end
        DistCI(i,:)=[Hrange(lg-max(lowind)+1)  Hrange(max(highind)+lg+1)]; % computing CI for distorted level
    end
end
end

%% Shock ordering constraint
function [c,ceq] = ordering(thetahat) % imposes ordering constraint: largest proportional variance change is in policy shock, epsilon 2
c = thetahat(4)/thetahat(2)-thetahat(5)/thetahat(3);
ceq = [];
end
